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Ando’s model provides a rigorous quantum-mechanical framework for electron-surface 
roughness scattering, based on the detailed roughness structure. We apply this 
method to metallic nanowires and improve the model introducing surface rough¬ 
ness distribution functions on a finite domain with analytical expressions for the 
average surface roughness matrix elements. This approach is valid for any rough¬ 
ness size and extends beyond the commonly used Prange-Nee approximation. The 
resistivity scaling is obtained from the self-consistent relaxation time solution of the 
Boltzmann transport equation and is compared to Prange-Nee’s approach and other 
known methods. The results show that a substantial drop in resistivity can be ob¬ 
tained for certain diameters by achieving a large momentum gap between Fermi level 
states with positive and negative momentum in the transport direction. 
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I. INTRODUCTION 


Modeling surface roughness (SR) scattering of electrons in metallic nanowires is quite 
challenging and often requires drastic approximations. Early attempts by FuchsP and later 
SondheimerP introduced a specularity parameter p (also known as the Fuchs parameter) 
that represents a probability for diffusive scattering at the boundary with the electron los¬ 
ing forward momentum. While predicting the experimentally observed 1/width behavior of 
the resistivity in thin metallic Elms or wiresp^this assumption artihcially imposes a bound¬ 
ary condition on the Boltzmann distribution function, neglecting the quantum-mechanical 
nature of the scattering events, the impact of conhnement and the detailed roughness pro- 
hle. It is therefore expected that the Fuchs-Sondheimer model, or the extension by Mayadas 
and Shatzkes to include grain boundary scattering,^ fails when the diameter of the metallic 
wire is reduced and becomes comparable with the roughness size or enters the regime where 
quantum-mechanical effects such as energy quantization become important. 

A predictive model for SR scattering is very important, because it is commonly used to 
extract information on the relative contribution to the overall resistivity, also coming from 
other scattering mechanisms such as phonons or grain boundaries. It is commonly stated 
that SR scattering is one of the dominant scattering mechanisms in metallic nanowires 
and therefore an important bottleneck for further downscaling of semiconductor devices, 
where they serve as interconnects. One should be careful distinguishing between resistivity 
contributions relying on the Fuchs-Sondheimer or Mayadas-Shatzkes model because they 
could be invalid for the dimensions of the wires in the experimental setup. The relative 
contribution of SR scattering to the resistivity could be properly confirmed if the results of 
simulations and experimental data would appear to be in agreement for the resistivity values 
and scaling with width (or diameter), without the need of non-physical fitting parameters 
such as the Fuchs parameter. A model that only depends on the details of SR (or its average 
properties) is therefore needed. 

Being a well-known approach developed to this end, Ando’s model employs Fermi’s golden 
rule to calculate the SR scattering probability, thereby relying merely on two statistical 
properties, the SR standard deviation and the SR correlation length.^. A commonly used 
approximation is the Prange-Nee approximation,^!^ replacing the overlap of wave functions 
in the conhnement directions by the value of the wave functions at the boundary times the 
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barrier height in the limit of an inhnite barrier. Recently, the resistivity was obtained for 
narrow metallic nanowires (width and height ranging from 1 to 6 nm), using the Prange- 
Nee approximation in Ando’s model. It was shown that the resistivity due to SR can have 
very different values and scaling behavior for different values of SR standard deviation and 
correlation length.^ 

A disadvantage of the Prange-Nee approximation is that it takes the inhnite barrier limit 
of an expansion for small roughness sizes, which often falls short to determine the wave 
function overlap, especially in the case of large roughness sizes, highly oscillating wave func¬ 
tions or low barrier heights. For metallic nanowires one should dehnitely be careful, as there 
are many subbands and, hence, high wave vector values in the conhnement direction are 
generally present. Quite recently, a new roughness model was proposed to improve the com¬ 
putation of the wave function overlap in thin-body fully depleted SOI device simulations, 
using a bivariate normal distribution function for the roughness statistics.!^ This method 
requires numerical integration for the large amount of subbands in metallic nanowires dras¬ 
tically slows down the simulation. Therefore we propose two variants of a novel approach, 
based on distribution functions on a hnite domain. In the regime of wire dimensions that 
can be simulated with all methods mentioned above, we will make the comparison. 

The paper is structured as follows. In section [IT] the SR formulation and statistics are 
briefly summarized. In section [ITT| this formulation is used to retrieve the SR matrix elements, 
adopting the Prange-Nee approximation, the bivariate normal distribution function and the 
two variants of the new approach, based on hnite domain distribution functions. We also 
briehy discuss how the SR matrix elements are used in the Boltzmann equation and the 
conductivity calculation, but this can be found in more detail in Moors et ah®! Results are 


given and discussed in section [IV] for a single subband toy model and more realistically sized 
metallic wires, followed by a conclusion in section jVj 


II. SURFACE ROUGHNESS 

Consider an ideal wire with length and rectangular cross section (width and height 
Ly), represented by a rectangular hnite potential well with barrier height U. The ideal wire 
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Hamiltonian Hq is then given by: 


IT ( \ - j if 0 ^ ^/y ^ 1^1 ^ 

^0 Hi - ^ .2 o ’ 

l-fev^ + c/, else 

with an effective mass rrie assumed to be appropriate for the conduction electrons in the 
metallic nanowire under consideration. Periodic boundary conditions in the transport di¬ 
rection 2 : are also assumed, neglecting any effect of the contacts. SR is modeled by SR 
functions that describe the deviation from the smooth, ideal boundary surface for each of 
the four boundaries: S^^^{y,z), z), S^^^{x,z) and {x, z). Typically one sup¬ 

poses Gaussian or exponential autocorrelation functions for the SR functions and in this 
paper we consider them to be Gaussian, 

{S^=\y, z)S^=\y, z)) = ( 1 ) 

with standard deviation A and correlation length A. We only consider SR around the 
ideal wire for each boundary separately thus neglecting edge effects of the autocorrelation 
functions at the intersection of two boundaries. The autocorrelation functions will be used 
in the following section to calculate all the SR matrix elements. 

III. TRANSITION PROBABILITY MATRIX ELEMENTS 

The probability for an electron to scatter from an initial state | i) to a hnal state | /) 
due to a perturbation Hamiltonian V is given by 

P(l *> ^1 /» = Y K* I r I - B /), (2) 

according to Fermi’s golden rule with Ej and Ef the energy of initial and hnal state re¬ 
spectively. The scattering probabilities enter the collision term of the Boltzmann transport 
equation (BTE) through the relaxation time approximation, such that the following equa¬ 
tions have to be solved to obtain the relaxation times: 

^ ^ nieLz I I f)f A _ /o'. 

where we assumed running wave solutions along the transport direction with wave vectors kf 
and for the initial and hnal state. The relaxation times can be calculated self-consistently. 
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... HS°^- [x-S--^(y,z),y,z] 
. V^=°{x,y,z) _ 


FIG. 1. The potential energy part of the ideal Hamiltonian, Hg°*'(r), is shown (continuous line) 
as a function of x for a fixed y and z, together with the Hamiltonian shifted by z) (dashed 

line) and their difference V^^^{x,y,z) (dotted line). To obtain the potential for SR of the x = 0 
boundary, only the part of close to x = 0 should be considered. This part is shown hatched 
on the figure. 

or by applying commonly used approximations for the ratio of relaxation times on the RHS. 
We only consider states at the Fermi level because SR scattering is considered to be elastic, 
while room temperature is very small compared to the Fermi level in metals,^ hence 
the deviation from the step function in the T = 0 Fermi-Dirac distribution is negligible. 
Using the nanowire model introduced in section [TT| a matrix element is in general given by: 

{i\V \ f) = ^ Jdrip* {x, [H{r) - Ho{r)] 'ipfix, (4) 

with 'ip{x,y) = 'ip{x)'ip{y) and 'ip{x/y) the finite potential well wave functions (the argument 
indicates the corresponding transversal direction). For SR at the x = 0 surface the Hamil¬ 
tonian H can be written as the unperturbed Hamiltonian shifted by the appropriate SR 
function: 


/C-" (x, y, z) = H„[x- z), y, z] . 

In this way the perturbation = ffx=o _ easily be obtained, as can be seen in 

Fig.lD The shift also acts around x = L^, which should not influence the x = 0 SR matrix 
element. Therefore the potential is always considered to be zero when |x| > \S^^^{y, z)\ 
while \S^^^{y, z)\ is assumed to be considerably smaller than for all y and 2 ; to exclude 
the unphysical situation in which the rough surface crosses the opposite wire boundary. The 
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expressions for the SR matrix elements coming from the other wire boundaries are similarly 
obtained. 


A. Prange-Nee approximation 


Instead of solving the integral in Eq. the Prange-Nee approximation expands the 
integrand up to hrst order of the SR function S^^^{y,z) while taking the inhnite barrier 
limit of the hrst order term.^ For the x = 0 surface, the expansion up to hrst order of the 
SR function yields: 


R' 


a :=0 


/) 




y 

'd!/ f dz z),p; (y) ,p, (y) (5) 

-L,I2 


The value of U'ip*{x)'ipf{x) at x = 0 depends on the height of the potential well and the 
corresponding values of the wave function at the boundary. The limit for U —?• -|-oo gives 
2^EfEJ/Lx with the energy of the initial/hnal state coming from the x-direction, 
which can also be written as {h?/2me) {d/dx)'ip*{x){d/dx)'ipf{x)\^^Q. To obtain an expression 
for the scattering probability Eq. the result of Eq. should be multiplied by its complex 
conjugate, making it possible to insert the correlation function of Eq. and calculating the 
scattering probability analytically® 

_ _ _ -\-Lzl‘2^ Ly -\-Lz/‘2. 

dy Jdz ( 6 ) 

0 - L ^/2 0 - L ^/2 




EfE^f A" 
[——- — 

LI LI 


dy / dz 


Assuming inhnite barrier height around the wire, the wave functions of the conhned direc¬ 
tions are given by \J2/L^ sin {jIx'kx/L^) with positive integer subband index (analog for 
y-direction). The length-dependent part in the limit S> A gives: 



(fcf-fcJ)2A2/8 


( 7 ) 
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while the remaining integration over y,y gives rise to an additional form factor such 
that: 


(|(i|W»|/)f)«4y| 


''2 LI L 


F^f = j^y j'^Ky)^fiyWfiy)My)e 

0 0 


The scattering probability is highly suppressed when \kf — kj\ > 2a/2/A. With this obser¬ 
vation we can define a critical momentum gap: 


A*:?,,.. <P», = 2^/A. 


( 8 ) 


The diminishing scattering probability for large Ak^ is a crucial observation also for the 
improved solution of the SR matrix elements. 


B. Improved solution 


The shortcomings of the Prange-Nee approximation become clear when we write the SR 
matrix element differently: 

Ly +L./2 S^=°{y,z) 


{ilV^ ^ \ f) = ^ jdy jdz xp*{y)'4)f{y)e j 

0 -L^/2 0 

The resulting expression in Eq. [^contains the integral (surface integral SI): 

S^=°{y,z) 

U [dx r^ (x) (x) ^ Shj [S^=\y, z)] , 


(9) 


( 10 ) 


which is a functional of the SR function. The integral is approximated in Eq. |^by the width 
of the integration region multiplied by the value of the integrand at x = 0. It is a good 
approximation as long as the integrand is not changing much over the integration region, 
namely if the SR deviation size is small enough. However, in the limit U —>■ -1-cxd the lowest 
order gives a finite value while the surface integral is divergent. We can see this more clearly 


by expanding Eq. 10 up to second order of the SR function: 


SIfj” « US--°{y, z) + U [S*-"!;/, z)f - 


7 




The same result is obtained by rewriting if with Hq in Eq. expanding it up to second 
order of the SR function and using partial integration to tackle the derivative acting on 
the Dirac delta function. If we now take the limit U —)■ +oo, the hrst term reaches a 
hnite limit but the second term diverges. It reflects the divergence of an integral over a 
hnite region with an integrand that is inhnitely large. The divergence was put under the 
rug by approximating the wave function by its value at the boundary, where the wave 
function should become zero when the potential is taken to be inhnite. This analysis shows 
that approximating the surface integral up to hrst order of the SR function is not fully 
compatible with the inhnite barrier limit. Therefore the validity of simulation results based 
on Prange-Nee can rightfully be questioned in some cases. Expanding Shj up to hrst order 
of the SR function without an inhnite barrier limit, thus requiring an extra parameter t/. 


might be a considerable improvement and is considered in section IV, referred to as the hrst 


order approximation. One could also expand Eq. 10 up to higher orders of the SR function 
to obtain a even better approximation. For realistic sizes of SR however, we expect this 
approach to be impractical because the wave functions can be highly oscillating (see Fig. 
and we would have to go to very high orders to capture the oscillations of the integrand. 

Instead of considering higher order expansions of the surface integral Shj, the integral 


in Eq. IT can be calculated analytically, taking into account the contribution of every order 
of the SR function. There are two diherent solutions, depending on the sign of the SR 
function, one being the solution for the product of exponential tails, the other for the 
product of oscillating wave functions inside the rectangular cross section of the wire. Below 
is an example for the surface integral at the a; = 0 snrface with initial state | i) and hnal 
state I /). We write the wave functions (only considering even wave functions with respect 
to X = La;/2 in the text, odd wave functions can be treated analogously) as: 


'ipi{x) = < 


for X < 0 

Bi cos [ki (x — La./2)], for 0 < X < La; • 
A.e-ai(x-L^)^ for X > La; 


The exponential tail on the other side of the wire in the integral of Eq. can safely be 
neglected because we assumed that the SR function is never larger than the wire width or 
height. The real constants Ai, Rj, ai and k^ are obtained from the hnite potential well 


solution in the x-direction. Insertion into Eq. 10 leads to: 









( 11 ) 


SI,J [S^=^] = I 

SI+^ ^ Y. 


SI+ [S^=°], > 0 


SI-/ 


, Sl-j[S^Y^-UAAf 


\ _ g(ai+a/)5'“ 


^x=0 ^ Q ' ‘’J ^ ^ ■ ■' Ui + ttf 

sin [{kf ± ki) L^/2] - sin [{kf ± ki) (L^ - 2S^=°) /2] 


^ kf±ki 

Next, these solutions have to be inserted into the SR matrix element given by Eq. [^and the 
average has to be taken over the absolute value squared of the matrix element appearing in 
the BTE (Eq. |^. We obtain for the contribution of the surface near x = 0: 


|(, I V—» I /)p) 


( 12 ) 


0 -L^/2 0 -L,/2 


d~z 




Although the SR of different boundary surfaces is considered to be uncorrelated, we do 
observe couplings between them in the average of the absolute squared of the matrix element. 
This is because SI [A] is in general not symmetric around S' = 0 thus leaving the average 
different from zero. Below the general form of the coupling between the x = 0 and y = 0 
surface is given: 


{{t\v^=^\f){f\vy=^\t)) (13) 

Ly -\-Lz}‘2 Lx ~\'Lzl‘2 

= J2 j^y (sI/7°SI^,7°) 

^ 0 -L^/2 0 -L^/2 


Inserting the expressions from Eq. in Eq. or Eq. the average over the SR 
functions can be split up in different regions depending on their sign: 

(14) 


Sh,/Sh,/) = (SI+^Sl7^) + 25? ((SI+^SI-^) j + (SI-^SI-^) , 
(SI./) = (sit/> + (SI-,). 


To continue, the average over the SR function in the surface integral should be calculated. 
This requires a distribution function f{S,S) that is consistent with the Gaussian autocorre¬ 
lation function of Eq. There are several ways to obtain an average that is consistent with 
the Gaussian autocorrelation and three approaches are discussed in the following subsections. 
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FIG. 2. The product of wave functions with the highest and next to highest k value is shown {Lx ~ 
2.15 nm) as the green curve, with the integral of Eq. 10 from x = 0 up to (a) 5^=0 = —0.5 nm (b) 


^a;=o _ being the green shaded area under the curve. The Prange-Nee approximation replaces 
the wave function product by the red dashed line, whereas the first order approximation replaces 
it by the blue dotted line, such that the area under the wave function product is approximated 
respectively by the green or blue shaded rectangular area. 


Bivariate normal model 


We can suppose a bivariate normal distribution function, as was done recently by Lizzit 
et al.^ For the x = 0 surface this gives: 

^2 + ^2 _ 2^-(r, f) 

2A2[1-G(r,f)'] J’ 

G(r,f) = 




= dS dS 


Sl-f [-5] Sl^f[S] 


exp 


27rA2A/l — C fr. rV 


As the latter integral cannot be evaluated analytically, it needs to be computed numerically 

over 


for every pair i, / and every possible value of 0 < C{r, r) < 1. The integration of Eq. 15 


the quadrant in the (S', S')-plane with S', S' > 0, denoted by the functional r)], can 

be htted as a power series of C'(r, f): 


C+f [C (r, f)] ^ 4’+ + (r, f) + (r, f) + ..., (16) 

with htting parameters c ^'12 ■ Fitting parameters can be found for the integration 

over every quadrant of the (S', S')-plane in the same way, with the superscripts ± and ± 
representing the sign of S and S respectively. The corresponding htting function can be 
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• Lx^O.ISnm ■ Lx=0.22nm 
^ Lx==0.25nm ^ Lx=*0.29nm 




5 

+ •-•■ 
U - 


(a) S' > 0, 5 > 0 


(b) S > 0, S < 0 



(c) S < 0, S > 0 


(d) S < 0, S < 0 


FIG. 3. Using Eq. 16 to fit the function Cjj(C') in the different quadrants of the (S, S)-plane 
for wire simulations with Lx ranging from 0.18 to 0.29 nm and A ~ 0.036 nm, A ~ 0.11 nm. 
These nanowire dimensions are unrealistic and relate to the single subband toy model discussed in 
section IIV A1 


plugged in the remaining integrals of e.g. Eq. 12, such that the remaining integration over 


the boundary coordinates can be done analytically. The htting function up to second order 
of C appears to capture the behavior quite well (see Fig. [^. 


For some of the pairs i, f, the integrand in Eq. [^is highly oscillating and the number of 
pairs increases rapidly for larger diameters, so it proves to be a tedious numerical job. We 
have not been able to optimize this method for metallic wire simulations with width or height 
up to 10 nm (see Fig. [^for a comparison of computation times). To overcome this problem 
of numerical integration, we propose two other distribution functions that are consistent 
with the Gaussian autocorrelation (for S and S SR functions of the same boundary) below. 
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FIG. 4. The finite domain for S and S is [— \/3A,+\/3A] x [— \/3A,+\/3A] and the weight of 
f{S, S) is given by 1 + 4(7 (r, f) /3 and 1 — 4(7 (r, f) /3 in the equal sign and unequal sign quadrants 
respectively. It is not difficult to check that this is consistent with the Gaussian autocorrelation 
function C(r,f). For S and S SR functions corresponding to different boundaries, C(r,f) is 
considered to be zero (no correlation). 

2. Finite domain model I 


An approach that allows for an analytical expression of the average over SR functions 
consists of considering uniformly distributed variables over a finite interval around the av¬ 
erage position of the wire boundary surface. The finite interval also makes sure that no 
outliers of the SR function are passing the other side of the wire, consistent with the as¬ 


sumption in Eq. 11 that a rough boundary does not cross the other side of the wire, unlike 


the bivariate normal distribution function in Eq. 15 A standard deviation and correlation 
length is introduced in the distribution function f{S,S) by changing the domain size and 
assigning different weights to different quadrants of the (S', S) plane respectively, as shown 
in Fig. 1^ while retaining the appropriate mean value of the SR functions ((S') = 0). Using 


this model for the average in Eq. 14, we get: 


SI.jSIi,/ ) = [1 + 4C (r, f) /31 ( SI.,/SIi,^ )"'* + [1 - 4C (r, f) /S] (SI.,^SI,, 


off-diag. 


( 17 ) 


(SI../SI,,, 


diag. 


1 

12A2 


/ 

VSA 

p 

2 

0 



dS Sll, IS] 

+ 

/ dS SI-^ [S] 


V 

J 
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J 

-+3A 
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SI./SI,, 


ofF-diag. 


12A2 


/ V3A 0 

JdSSi:jl.S] JdS SI-; 

V 0 -V3A 




0 \/3A 

+ /dssr,|s| fdssit 


ij 


S 


-x/3A 
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When looking at Eq. 11, it is clear that only an integration over exponentials with complex 


argument remains, while applying the bivariate normal distribution (Eq. 15) leads to a 
remaining integral over a product of exponential, Gaussian and error functions for which an 
analytical expression could not be found. 

A downside of this approach is the unphysical negative weight that is attributed to the 
anti-correlated quadrants when ^(r, f) >3/4. This should be kept in mind when using this 
approach to simplify the integrals over the SR functions. 


3. Finite domain model II 


Instead of assigning different weights to different quadrants, the model proposed in this 
section consists of a sum of a fully correlated part and a fully uncorrelated part. This also 
makes the average over the SR functions analytically solvable: 


= C(r,f) (SlySI.,+ |1 - C(r,f)l (SI.,/SI.,^ 

n/3A 


SI./Shj 


vd2A 


d^|Shj[^]^ (Sh,/S\; 


12A2 


-V3A 


(18) 

x/3A ^ 

jdSShj[S] . 

-viiA 


The uncorrelated part has a uniform distribution in the {S, S')-plane, whereas the correlated 
part has a distribution function that only allows identical values of S and S, as shown in 
Fig. 1^ The full SR matrix element expression of this model and hnite domain model I can 
be found in appendix [A} 

Even though this model and the hnite domain model I are consistent with the bivariate 
normal distribution function for (S'™S'") with m -|- n < 2, the higher moments (m -|- n > 2) 
will be different. The higher moments of this model and the other models are compared in 
Fig. § They are different for each distribution function and the differences become more 
substantial when higher values of n or m are considered. The hnite domain models agree 
quite well but the expectation values deviate substantially from the bivariate normal model. 
The expectation values are in general larger for the latter because higher values of the SR 
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FIG. 5. The finite domain for S and S is [—\/3A,+\/3A] x [—\/3A,+\/3A] and the weight of 
this domain is 1 — C (r,f). The correlated part is the diagonal S = S with —-v/SA < S < ^/SA 
and weight C (r, f). This model is also consistent with the Gaussian autocorrelation function. For 
S and S SR functions corresponding to different boundaries, C (r, f) is considered to be zero. 

function contribute (infinite domain). All the distribution functions respect the symmetry 
and give zero when n + m is odd. 

IV. RESULTS & DISCUSSION 
A. Single subband toy model 

A test case model for the finite domain models proposed in previous section is a wire with 
only a single subband crossing the Fermi level. This can be realized with a low effective mass, 
conduction electron density or diameter. While in principle being possible (and already being 
proposed and investigatecP®^!!) for semiconductor nanowires, it is practically impossible to 
realize for metal nanowires. It serves as a good test case however because there are only two 
states at the Fermi level hence scattering is limited to Umklapp scattering and we can clearly 
see the differences of the different methods and the effect of the momentum gap Ak^ on the 
transport. We simulate the single subband model by considering a very narrow wire and 
by fixing the Fermi level with the bulk electron density. In Fig. [^the resistivity is shown, 
using all the methods introduced in section [TTj Fig. [^clearly shows the large overestimation 
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(a) 0 < C < 1 


(b) C = 1 



(c) C = 0.25 


(d) C = 0.75 


FIG. 6. The expectation values of products of SR functions are shown for several values of corre¬ 


lation C, using the different distribution functions that are proposed in section III B 


of the Prange-Nee model, more or less by an order of magnitude. The results are very 
different because we are simulating an extremely narrow nanowire. The inhnite barrier 
approximation is not working well because there is substantial wave function penetration 
into the potential barrier region. The other methods are much better in agreement, with very 
good agreement between the bivariate normal model and hnite domain model I. The hrst 
order approximation of the wave function surface integral does also not differ much different 
from the former two models. The hrst order approximation is expected to work well for this 
toy model because there is only one subband and therefore the wave function, having a large 
wave vector, varies little over the SR standard deviation. A hrst order approximation of the 
integral in Eq. |^is acceptable in this case. The approximation is expected to be worse when 
more subbands (and higher wave vectors) come into play. 

The increase of resistivity p for decreasing Ak^ seems to follow a p = pmax. exp (—LAk^) 
scaling quite well for the Prange-Nee model, but the other methods deviate from this with an 
even larger suppression for large Ak^. The parameter pmax. dihers from method to method. 
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^PN 

^ First order 
^ Finite domain I 
Finite domain II 
^ Bivariate normal 


FIG. 7. In (a) the dispersion relation of a single subband toy model is shown with the momentum 
gap Ak^ indicated. The resistivity values obtained with SR matrix elements are shown in (b) 
in a log scale, using the Prange-Nee model (red triangles), first order model (purple upside-down 
triangles), finite domain models (I: blue circles, II: yellow squares) and finally the bivariate normal 
model (green diamonds). The simulated wire has a square cross section with width and height 
ranging from approximately 0.25 nm to 0.29 nm and roughness properties are: A ~ 0.036 nm, 
A ~ 0.11 nm. The Fermi level is obtained by fixing the density to the bulk conduction electron 
density in copper (7 eV in bulk), but together with the very low width and height, fixing the Fermi 
level in this way is just a tool to obtain a single subband model and does not represent a realistic 
nanowire. 

but all the methods agree quite well on a parameter L 0.1 nm. Although the scaling of 
Eq. is different from the scaling in the simulation results, both indicate an exponential 
suppression for large Ak^ and it gives a critical length scale 0.04 nm of the 

same order of magnitude as L obtained from the £t. This behavior is very important because 
every wire is characterized by a minimal Ak^ between positive and negative wave vector 
states and it appears to have a crucial impact on the self-consistent solution of relaxation 
times when more subbands are considered. Hence, it also crucially affects the resistivity, as 
will be discussed in more detail in subsection IIVBI 

Other than the resistivity values of the different methods, there is also the difference in 
computational burden. We expect the Prange-Nee model to be the fastest and the numerical 
integration of the bivariate normal model to be the slowest. This is confirmed by comparing 
the simulation times in Fig. Currently the simulations of the hnite domain models are 
more or less one order of magnitude slower than the Prange-Nee model, while the bivariate 
normal model performs even slower by two orders of magnitude. While it must be said that 
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FIG. 8. The simulation time for the different methods used to calculate the SR matrix elements is 
shown on a logarithmic scale (in milliseconds) for the single subband simulation of the four data 
points presented in Fig. (b). The simulation is performed on a single i5 3.4GHz CPU and is 
based on a global adaptive numerical simulation strategy. A realistic nanowire size would increase 
the number of computations by a factor of lO'^ to 10®. 

our implementation of the numerical integration might be further optimized to enhance the 
computation speed of the latter, it will be hard to beat the performance of simulations using 
hnite domain model I or II because they only consist of evaluations of analytical expressions. 

B. Metallic nanowires 

In order to simulate realistic metallic nanowires, the model parameters have to be opti¬ 
mally determined. For the barrier height we consider the value corresponding to vacuum 
using the relation U = + W, with W the work function and Bp the Fermi energy for 

bulk.l^ The states of the subbands that cross the Fermi level are obtained self-consistently 
for every wire dimension by hlling up all the subbands up to the appropriate carrier den¬ 
sity of the considered material. All the simulations rely on the underlying assumption of a 
close to rectangular cross section and an (isotropic) effective mass for the conduction elec¬ 
trons (taken from measured bulk values), which can rightfully be questioned for very narrow 
metallic nanowires and some materials. However, the presented approach with hnite domain 
distribution functions and their analytical expressions are still applicable to more detailed 
wave functions and subband structures (obtained with ab initio methods for example), as 
long as the wire boundary consists of hat surfaces. In this case, one can dehne SR functions 
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on those surfaces and apply Eq. m straightforwardly by expanding the wave functions 
on a basis of exponentials or sines and cosines.^ 

The resistivities, obtained with the different methods and the self-consistent solution of 
the BTE, are compared in Fig. In contrast to the one-band toy model, hnite domain 
model I now differs most from the other methods and Prange-Nee works better, although 
it does not follow the same trend as all other methods. The first order approximation and 
hnite domain model II now seem to be closest in agreement, with only small difference in 
resistivity value for all diameters. The overall trend for decreasing diameters is an increase 
of the resistivity (not with the Prange-Nee model results), but the results show a very low 
resistivity value for the wire with sides between 3 and 3.5 nm for all methods. Comparing 
the very low resistivity values obtained with the different methods, we observe quite big 
differences. This can be expected because the important states for back-scattering have 
quite high wave vectors in the conhnement directions and the different methods compute 
the wave function overlap due to roughness quite differently for these states. The resistivity 
drop appears to be related to the size of the minimal momentum gap between the positive and 
negative state with lowest absolute value. When this gap exceeds the critical momentum 
gap (Eq.g, back-scattering is highly suppressed and only scattering to other states without 


substantial loss of forward momentum or current takes place. This can be seen in Fig. 11 


where the band structure and relaxation time of the two lowest diameter data points of 


Fig. 10 are shown. The momentum gap of the larger diameter is approximately , 

which increases the lifetime of the low states drastically. 

Another interesting aspect of the resistivity drop is that it is only visible in the self- 
consistent solution of Eq. This solution is shown in Fig. [To] together with the approximated 
solution by putting the relaxation time ratio on the right-hand-side equal to 1, assuming 
equal lifetimes for the initial and hnal states.!^ The resistivity drop is only apparent in the 
self-consistent solution of the BTE because it is the only solution that solves the coupled 
system of equations, determining the scattering lifetimes correctly. Two data points of 


Fig. I^are analyzed in more detail in Fig. IT From that figure we can estimate the critical 
gap and it appears to agree with 2 nm“^, obtained from Eq. When 

the incoming scattering is neglected in Eq. by putting the relaxation time ratio equal to 
zero, the relaxation times are largely underestimated and the predicted resistivity is largely 
overestimated. 
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FIG. 9. The resistivity for copper nanowires with square cross section of sides ranging from 2 
to 7 nm (D = Lx = Ly, A = 2acu ~ 0.7 nm, A = 5acu ~ 1.75 nm) are shown here, using the 
different methods explained in section]^ Prange-Nee model (blue circles), first order model from 
Eq. (yellow squares), finite domain model I (green diamonds) and finite domain model II (red 
triangles). 



FIG. 10. The resistivity for copper nanowires with square cross section of sides ranging from 
2 to 7 nm {D = Lx = Ly) are shown here, using the self-consistent and approximate solution, 
substituting the relaxation time ratio on the RHS of Eq. by 1, of the BTE with SR matrix 
elements obtained from finite domain model I (A = 2acu ~ 0.7 nm, A = 5acu ~ 1.75 nm). 

Finally, suppose that the momentum gap is around 1 or 2 nm“^, typically the case for 
nanowires with width and height up to 10 nm. Then the critical correlation length A^rit. 
will be of the order of 1 nm. If the SR correlation length is substantially larger than 
Acrit., the scattering matrix element will be exponentially suppressed by Eq. (or a similar 
expression for hnite barrier heights) and the resistivity contribution becomes very small, 
even if the SR standard deviation is quite large. This is in agreement with results reported 
by Leunissen et al.^ who found that the influence of roughness added to a wire has no 
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/t^(1/nm) 


(a) Subbands D « 2.15 nm 


(b) Subbands D « 3.25 nm 



/d{1/nm) ;d(1/nm) 


(c) Relaxation times D « 2.15 nm 


(d) Relaxation times D « 3.25 nm 


FIG. 11. All the subbands that cross the Fermi level of a copper wire with Lx = Ly equal to (a) 
6 Ocu ~ 2.15 nm, (b) 9 Ocu ~ 3.25 nm are shown as a function of . The self-consistent (green 
circles) and approximate (relaxation time ratio = 0: yellow squares, relaxation time ratio = 1: 
purple diamonds) relaxation time solutions of Eq. with SR matrix elements, obtained with finite 
domain model I (A = 2acu ~ 0.7 nm, A = 5acu ~ 1-75 nm), are shown for every Fermi level state 
with positive k^ in femtosecond log scale in (c-d). The red dashed line is showing the Fermi wave 
vector kp = \/2meE'p/fi. 


impact. The critical correlation length is on an atomic scale, rather than the scale of the 
wire dimensions or the scale related to the etch process (typically larger than the atomic 
scale). Ronghness properties on atomic scale have already been considered in semicondnctor 
device modeling, experimentally verihed for those types of bonndaries with AFM 
and TEM.I2SH2I1. The SR on atomic scale appears to be crncial for metallic nanowire modeling 
as well and it is very nnlikely to be much dependent on wire width or height. Therefore the 
SR standard deviation and correlation length were considered to be independent of the wire 
dimensions in the simulations presented above. 
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V. CONCLUSION 


The investigation of SR scattering, partially responsible for an increased resistivity in 
narrow metallic nanowires, requires a good model for comparison with experimental data. 
We propose a model based on Ando’s model to model SR scattering exploiting the statistical 
features of roughness, rather than more empirical approaches like the Fuchs-Sondheimer 
model. With this model we can compare experimental data for metallic nanowires both 
on the level of resistivity measurements and detailed SR properties. This allows for proper 
testing of the SR model as well as improved resistivity predictions and simulations. 

A crucial observation is that the BTE under the relaxation time approximation should be 
solved self-consistently to obtain the relaxation times and resulting resistivity correctly for 
SR scattering. While this method was already proposed in Moors et ah, ^ the Prange-Nee 
approximation was always used before to obtain the SR matrix elements that are based 
on Ando’s model. We have now shown that the Prange-Nee approximation has serious 
shortcomings by neglecting the wave function oscillations and considering an infinite barrier 
limit. Because the wave functions can be highly oscillating, an alternative approach has 
been presented here under two variants to obtain Ando’s matrix elements. The two variants 
are consistent with Gaussian SR autocorrelation functions and allow for a fast and accurate 
evaluation of the SR matrix elements with any roughness size. The statistical properties of 
surface roughness are captured in both cases by an appropriate distribution function on a 
hnite domain. This new approach has also been compared with a recent proposal that also 
gives the SR matrix element without assuming small roughness sizes, but requires numerical 
integration when a bivariate normal distribution function is invoked. 

Comparing all the different methods, the Prange-Nee model surprisingly works quite well 
for resistivity estimation of metallic nanowires with sides of a few nanometers. It clearly gives 
wrong results in the single subband toy model due to a large difference between the finite 
and infinite potential well solutions. In general the first order model works much better, even 
though the oscillations of the wave functions normal to the boundary plane are completely 
neglected. The newly introduced methods are in good agreement for the resistivity scaling 
trend and are faster than the method based on a bivariate normal distribution function by 
two orders of magnitude. 

An important observation in this work is a resistivity drop for certain diameters, con- 
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firmed by all the self-consistent BTE solntion with all the different methods. The crncial 
parameter is the wave vector gap between the Fermi level states with negative and positive 
wave nnmbers. An exponential snppression of cnrrent-loss indncing scattering occurs when 
this gap exceeds a critical inverse length scale that is found to be proportional to the inverse 
of the SR correlation length for all the simulations. This allows for a conductivity enhance¬ 
ment by engineering the subbands in such a way that this gap is maximized, a proposal that 
was already made for semiconductor nanowires by Sakaki in the SO’s.^^^® 
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Appendix A: Distribution functions on finite domain 


The expression for the average SR matrix element, using the distribution functions pro¬ 


posed in section 111 B 2 and section 111 B 3, is quite involved and it is easy to loose track 


of the calculation. The full expressions for the matrix elements containing the four rough 
boundaries of a wire (with equal distribution functions and mutually uncorrelated) that were 
used in the simulations are given below for reference: 
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Note that even though the integration results are left implicit, every integral in the expres¬ 
sions can be written as an analytical expression. 


25 



